Molecular Dynamics Study of Charged Dendrimers in Salt-Free Solution: 

Effect of Counterions 
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Polyamidoamine (PAMAM) dendrimers, being protonated under physiological conditions, repre- 
sent a promising class of nonviral, nano-sized vectors for drug and gene delivery. We performed 
extensive molecular dynamics simulations of a generic model dendrimer in a salt-free solution with 
dendrimer's terminal beads positively charged. Solvent molecules as well as counterions were ex- 
plicitly included as interacting beads. We find that the size of the charged dendrimer depends 
non-monotonically on the strength of electrostatic interactions demonstrating a maximum when the 
Bjerrum length equals the diameter of a bead. Many other structural and dynamic characteristics 
of charged dendrimers are also found to follow this pattern. We address such a behavior to the 
interplay between repulsive interactions of the charged terminal beads and their attractive inter- 
actions with oppositely charged counterions. The former favors swelling at small Bjerrum lengths 
and the latter promotes counterion condensation. Thus, counterions can have a dramatic effect on 
the structure and dynamics of charged dendrimers and, under certain conditions, cannot be treated 
implicitly. 



I. INTRODUCTION 



Dendrimers, characterized by regular branching with 
radial symmetry, represent a unique class of poly- 
mers possessing a large number of functional terminal 
groups in their outermost dendritic shelli These "dream 
molecules'— have recently attracted a lot of attention due 
to potential applications in technology and medicine, see, 
for example, Refs. ETll lsl and references therein. Polyami- 
doamine (PAMAM) dendrimers, in particular, are com- 
monly used molecules in this field. They are proto- 
nated under physiological conditions?^ thus their high 
surface charge coupled to their ability to adapt con- 
formation in response to changes in their surrounding 
environment make PAMAM dendrimers very attractive 
candidates for biomedical applications such as drug and 
gene deliveryi&L&iLiSiiiiiSiiii Thus far, most theoretical 
and computational investigations have focused on neu- 
tral dendrimers^*iii 

Since the pioneering paper of Welch and 
Muthukumar^i there have been only a few compu- 



tational studies of charged model dendrimers ^LiLi2ii2i 
While these studies have provided a great deal of insight 
into the properties of dendrimers, it is obvious that their 
scope has been limited due to a major computational 
cost associated with modeling of complex dendrimers 
in a hydrodynamic solvent. Hence, the previous studies 
have treated the groups of a charged dendrimer using the 
linearized Poisson-Boltzmann (Debye-Hiickel) theory, 
and the solvent as a continuous medium. Despite these 
limitations, the above simulations demonstrated that 
the strength of electrostatic interactions may have a 
strong effect on the dendrimer's size as well as on other 
structural and dynamic properties. In particular, it 
was shown that the size of any charged dendrimer in- 
creases for an increasing effective screening length which 
corresponds to a decreasing salt concentration li&iLiSiiS 
As for more accurate treatments, Lee et aliSS simu- 
lated PAMAM dendrimers at an atomic resolution. They 
found that a PAMAM dendrimer swells significantly 
when pH decreases from a high pH value (no amines 
are protonated) to neutral (all the primary amines are 
protonated), and to low pH (all the primary amines and 
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tertiary amines are protonated) . However, they treated 
the counterions as well as the solvent implicitly (except 
for dendrimers of second generation which were solvated 
in explicit water) m& A more realistic atomic-scale sim- 
ulation of PAMAM dendrimers was performed very re- 
cently by Goddard et ai.j 2 !^ who accounted for water 
molecules and counterions in an explicit manner. They 
showed that the presence of solvent leads to a swelling 
of the dendrimer, while the effect of pH was found to 
be similar to that reported in earlier studiesi^S As a 
part of counterions were found to condense onto the 
dendrimerj2i it was stated that the presence of counte- 
rions increases the swelling of a dendrimer. The effect 
of counterions on dendrimer charges and related struc- 
tural aspects were not clarified, though, since for the 
counterion concentrations they used, the amine charges 
were strongly screened: the corresponding Debye length 
was about 3-4 A and hence close to the minimal value 
considered by Welch and Muthukumarii^ This implies 
that the effect of smaller counterion concentrations and, 
hence, the overall role of counterions has remained unre- 
solved. 

In general, there is reason to emphasize that the ap- 
plicability of the Debye-Hiickel approximation for den- 
drimers should not be taken for granted. It is valid in 
the dilute limit-here it means the situation where the 
separation between charges exceeds the Bjerrum length. 
However, that condition breaks down easily in dense 
structures such as dendrimers, which implies that in 
most of the cases counterions must be taken explicitly 
into account. Computational studies of flexible linear 
polyelectrolytes have demonstrated^ 3 * 2 ^^ that the pres- 
ence of explicit counterions is indeed crucial: increas- 
ing the strength of electrostatic interactions can provoke 
counterion condensation and, therefore, lead to dramatic 
changes in the polymer's structure and dynamics. It is 
reasonable to expect similar counterion effects also for 
dendrimers. A recent computational study of model den- 
drimers suggests that this is indeed the casei^i 

In this study we aim to clarify the role of counteri- 
ons. We have performed extensive molecular dynamics 
simulations using a generic cationic model dendrimer of 
the fourth generation. Counterions and solvent molecules 
are explicitly included. A systematic variation of the 
strength of electrostatic interactions have allowed us to 
examine the influence of counterions on the structural 
and dynamical properties of charged dendrimers. It turns 
out that the size of charged dendrimers depends non- 
monotonically on the strength of electrostatic interac- 
tions in the system clearly indicating a pronounced effect 
of counterion condensation. 



II. MODEL AND SIMULATION DETAILS 

We performed MD simulations of dendrimers of the 
4th generation using the freely jointed "bead-and-spring" 
model. Each bead of the last, or the terminal, generation 




FIG. 1: A sketch of the dendritic topology used in this study. 
Beads of the last, or the terminal, generation shell (the illus- 
tration is for a dendrimer of the second generation) are posi- 
tively charged (black). In the simulations we used dendrimers 
of the fourth generation. 

was set to carry a positive unit charge. The dendrimer 
was solvated in a box of non-polar molecules (beads) and 
an appropriate number of counterions was added to pre- 
serve overall charge neutrality. 

The bonds were modeled as harmonic springs with 
their potential energy given by U = ^k(l — lo) 2 , where 
lo is the reference bond length. The spring constant 
k = 10 3 x k^T /Iq was chosen to be rather high as 
has been done in earlier studies^ All short-range non- 
electrostatic non-bonded interactions were described by 
the Lennard-Jones potential [/ L j = 4e[(er/r) 12 — (a/r) 6 ] 
with e = O.SfceT. The bead-bead interaction diame- 
ter a — 0.8 lo was taken to be the same for all beads in 
the system. The Lennard-Jones potential was cut off at 
r c = 2 1 / 6 cr and shifted to zeroi^ 3 . The dendrimer was in 
a good solvent. 

The long-range electrostatic interactions between the 
charged terminal beads (with a positive unit charge 
on each) and negatively charged counterions were han- 
dled using the particle-mesh Ewald (PME) method; 28 i 29 
which has been shown to perform well in soft matter 
simulations The strength of electrostatic interac- 
tions is characterized by the so-called Bjerrum length 

AB = A T ' W 

47re e s fc B T 

where e s is the dielectric permittivity of the solvent. To 
examine the effects of counterions on the dendrimer's 
structure and dynamics, the Bjerrum length Ab was sys- 
tematically varied from to 8.0^o. This change in Ab 
may be interpreted as a variation of the dielectric prop- 
erties of the solvent (Ab ~ l/e s )- In ah, we studied 10 
different systems with Ab/?o taken to be 0.0; 0.2; 0.4; 0.6; 
0.8; 1.0; 2.0; 3.2; 6.4; and 8.0. 

The dendritic topology employed in our study is 
schematically depicted in Fig. ^ We studied a trifunc- 
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tional dendrimer of the fourth generation with single 
bonds between the branch points. The dendrimer had 
94 beads of which 48 were terminal ones. To generate 
initial configurations, the procedure of Murat and Grest 
was used^ The dendrimer was then solvated in a cu- 
bic box with 4692 solvent particles. That corresponds 
to a polymer fraction of about 0.02. The density was 
set to 1.688 In 3 (or, alternatively, to 0.864er -3 ), similar 
to Refs. 113 and |3 The linear size of the simulation box 
size was set to 14.16 Iq, and periodic boundary conditions 
were applied in all three dimensions. The ratio of the av- 
erage dendrimer size to the simulation box size was about 
0.2. 

For simulating charged dendrimers, unit positive 
charges were assigned to all beads of the terminal G4 
generation shell of the neutral dendrimer. This corre- 
sponds to physiological (neutral) pH conditions when 
the primary amines of a PAMAM dendrimer become 
protonatedpi To keep the system electroneutral, ran- 
domly chosen 48 solvent beads were converted to coun- 
terions by assigning negative unit charges on them. This 
system was used as the starting configuration for MD 
simulations of charged dendrimers. 

All systems with different Bjerrum lengths were first 
equilibrated for 10 6 time steps in the NVT ensemble with 
the time step At = 2.4 x 10~ 3 To, where t = l^m/e) 1 / 2 
is the characteristic time of the model (to is the mass of 
a bead). Temperature was controlled using the Berend- 
sen scheme^ while all the bond lengths were kept con- 
stant by using the SHAKE algorithm. 36 After equilibra- 
tion, data were collected from production runs of 10 
time steps each in the NVE ensemble. The time step in 
the production runs was set to 1.2 x 10 -3 tq (constraints 
were not applied to dendrimer 's bonds during the pro- 
duction). All the simulations were performed using the 
Gromacs package^2^ Each simulation was run on a sin- 
gle 3.2 GHz Pcntium-4 processor. 

To examine possible finite size effects, we performed 
additional simulation runs using simulations containing 
twice the number of solvent molecules. We examined 
the box size effects for a neutral dendrimer as well as 
for the charged dendrimers with Ab = 0.4 Z and Ab = 
6.4 Iq. For the case of neutral dendrimers we did not find 
any noticeable deviations when the larger box size was 
employed. For charged dendrimers, the average size of 
the counterion cloud around the dendrimer was found 
to naturally increase upon increasing the box size. This, 
however, did not produce any considerable changes in the 
structural or dynamic properties of dendrimers. 



III. STRUCTURAL PROPERTIES 

A. Equilibration 

As in previous computational studies^iSi initial equi- 
libration was monitored through the time evolution of 
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FIG. 2: Average mean-square radius of gyration of the 
whole dendrimer (solid circles) , and of the beads of the termi- 
nal dendrimer shell (open circles) as a function of the Bjerrum 
length Ab- Error bars are of the same size as the symbols. The 
error bars were estimated as the standard errors of mean by 
splitting trajectories into 10 pieces of 10 6 time steps each. 

the mean-square radius of gyration 
1 - 

R2 s=N E^-^m) 2 , (2) 

i=l 

where N is the number of beads and Rqm is the position 
of the center of mass (CM) of the dendrimer. This quan- 
tity, being a measure of the overall size of the dendrimer, 
was found to remain stable during the entire production 
run. 

The above, however, is not enough for a system con- 
taining unscreened charges. Is it further required that 
the ions have reached equilibrium. To monitor the equi- 
libration of ions, we first calculated the radial distribu- 
tion functions (RDFs) for ion-ion, ion-dendrimer bead, 
and ion-solvent molecule pairs (data not shown). The 
position of the first minimum defines the radius of the 
first coordination shell of a counterion. This procedure 
has been used in other related studies as wel L 39 i 40 ' 41 The 
radius was found to be ^1.3 Zo- The coordination num- 
bers were then calculated by counting the total numbers 
of beads within the first coordination shell. The coordi- 
nation numbers are found to become stable as a function 
of time during the equilibration period (see Sect. ITTj l. 

B. Dendrimer size and counterion condensation 

Next, we describe the structural characteristics of 
charged dendrimers at different strengths of electrostatic 
interactions. In Fig. we plot the average mean-square 
radius of gyration (Rg) of the dendrimer as a function of 
Bjerrum length Ab- Remarkably, the dependence of (RT) 
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FIG. 3: (Color online). Typical configurations of charged 
dendrimers at Ab = 0.8 lo (top) and Ab = 6.4 lo (bottom). 
Counterions are shown as spheres. Solvent is omitted for clar- 
ity. 



on Ab is found to be non-monotonic. Increasing the Bjer- 
rum length from zero leads to a pronounced swelling of 
the dendrimer, the maximum being at Ab — 0.8 lo as 
shown in Fig. [21 Further increase is found to have an op- 
posite effect, i.e., the dendrimer shrinks. Such a behavior 
resembles the findings for linearis and dendritic 26 poly- 
electrolytes, and is associated with counterion condensa- 
tion. Figure |3 shows snapshots from our simulations for 
the swollen and shrunk dendrimers. 

The As-dependence for the radius of gyration calcu- 
lated for charged beads of the terminal shell ((i?g)term) 
follows closely that of the whole dendrimer ((-Rg)dend)> 
see Fig. Interestingly, the (-Rg)tcrm exceeds the ra- 
dius of gyration of the dendrimer at all Bjcrrum lengths. 
For neutral dendrimers this finding is in agreement 
with experimental data^ and with Brownian dynamics 
simulations^ As demonstrated in Ref. I43L and as we 
proceed to show, the fact that (i?g)tcrm > (-Rg)dond does 
not mean that beads of the terminal generation shell are 
localized near the dendrimer's periphery. 

To clarify the role of counterions, we calculated the 




FIG. 4: Top: Average mean-square radius of gyration (Rg) 
calculated for counterions with respect to the center of mass 
(CM) of the dendrimer. Bottom: Average fraction of coun- 
terions condensed onto the charged dendrimer as a function 
of the Bjerrum length Ab- See text for the criterion used to 
identify condensed ions. 

mean-square radius of gyration (i?g)i on s of the counteri- 
ons, see Fig.0](top). That serves as a measure for the size 
of the counterion cloud around the dendrimer. It turns 
out that the ion cloud becomes smaller as the Bjerrum 
length increases. This decrease, however, is found to be 
monotonic. Furthermore, in Fig. 0] (bottom) we estimate 
the fraction of counterions condensed onto the dendrimer. 
For doing that we used a simple criterion, namely we 
counted counterions which had dendrimer beads in their 
first coordination shells. As Fig. 0] shows, the number of 
condensed ions increases continuously with the strength 
of electrostatic interactions in the system. 

The above suggests that there is a subtle interplay 
between two opposite effects upon increasing Ab: (*) A 
growth in the strength of the repulsive electrostatic in- 
teractions between the charged terminal beads and (ii) 
a strengthening of condensation of counterions due to 
the attractive interactions between counterions and the 
oppositely charged dendrimer beads. The former leads 
to a swelling of the charged dendrimer and the latter 
screens positive charges of the terminal dendrimer gener- 
ation giving rise to shrinking of the dendrimer. The max- 
imum of the dendrimer's size is found to be at Ab — 0.8 
(i.e. at Ab — cr), see Fig- EI This corresponds to swelling 
by ~ 12.5 % with respect to the neutral case. 
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dendrimer's periphery making the interior of the den- 
drimer significantly more accessible for solvent and coun- 
terions. 

Further increase in Bjerrum length leads to a shrinkage 
since counterion condensation starts to dominate. There- 
fore, the dendrimer's density profile changes towards a 
similar profile as observed at small Ab with one impor- 
tant exception: the amount of condensed ions increases 
considerably, and their density profiles overlap more and 
more with the density profiles of the charged terminal 
beads of the dendrimer (up to almost complete coinci- 
dence of the profiles in the case of extremely large Bjer- 
rum lengths, see Fig.|Hlfor Ab = 8.0 Iq). This shrinkage is 
accompanied by a reduction of voids in the dendrimer's 
interior leading to squeezing of solvent molecules out of 
the dendrimer, see Figs. El and [7| 



FIG. 5: Average coordination numbers of ions with the den- 
drimer's beads (solid circles), with the charged beads of the 
G4 shell (open circles), and with solvent beads (open squares) 
as a function of the Bjerrum length Ab. 

Condensation of counterions onto the dendrimer 
should be accompanied by a decrease in its hydration, 
which is characterized by the number of solvent molecules 
in the ion's first coordination shell. Indeed, that can 
be seen by computing the average coordination numbers 
of counterions with dendrimer's beads and with solvent 
molecules, see Fig. |SJ The increase of the strength of 
electrostatic interactions in the system provokes a con- 
densation of counterions and, correspondingly, leads to a 
loss of solvent molecules from the ions' first coordination 
(hydration) shell (Fig.JSJ. 



C. Number density profiles 

To analyze the locations of the dendrimer's beads, sol- 
vent molecules, and counterions in more detail, we cal- 
culated componentwise number density profiles p(r) as 
a function of the radial distance r from the center of 
mass (CM) of the dendrimer, see Fig. (Fig. [7| shows 
the corresponding radial bead distributions). For a neu- 
tral dendrimer the density profiles are in agreement with 
earlier molecular dynamics studies of dendrimers in an 
explicit solvent^ The overall number density of a den- 
drimer develops a plateau at r > 0.5 Iq from its CM. The 
beads of the terminal shell are broadly distributed, im- 
plying a considerable degree of back- folding. The solvent 
molecules penetrate deeply into the dendrimer and their 
density profile has a plateau at r > 0.5 Iq, see Fig.0 

At small Bjerrum lengths (Ab < 0.8 Zo) the charged 
dendrimer swells and the density profile plateau broadens 
considerably with a simultaneous drop in the plateau's 
magnitude, see, e.g., Fig. El for Ab = 0.8 lo. The charged 
beads of the terminal shell are re-distributed towards the 



D. Charge density profiles 

Next, we take a look at the screening of the dendrimer's 
charges by counterions. In Fig. |H1 we plot the charge 
density profiles p e (r) for the whole system as a func- 
tion of the radial distance r from the CM of the den- 
drimer. We recall that since our solvent is non-polar, the 
charge densities p e (r) are defined only by the interplay 
between the positively charged terminal beads and the 
negatively charged counterions. Figure |S] clearly demon- 
strates screening: the height of the positive /9 e (r)-peaks 
drops drastically with increasing Bjerrum length as com- 
pared to the /3 e (r)-peak for Ab = 0.2 Iq when the counte- 
rion condensation has just started. 

While the main peaks of all p e (r)-curves are located 
in the domain of the positive charges, starting with 
Ab = 0.6 Iq, one observes small minima close to the 
CM of the dendrimer where p e (r) becomes negative, see 
Fig. [8] We address these minima to the fact that some 
of condensed counterions can be localized deep in the 
dendrimer's interior rather than beside charged terminal 
beads giving rise to an excess of negative charges close to 
the dendrimer's center of mass, see Figs.Eland[7| In turn, 
on the right of p e (r)-peaks we see another region where 
the overall charge of the system is also negative. This re- 
gion of negative charge, being associated with "unbound" 
counterions in bulk solvent, becomes smaller with Ab and 
eventually disappears at large Bjerrum lengths. 



IV. DYNAMIC PROPERTIES 

The molecular dynamics simulation technique provides 
an access to the dynamic properties, and next we discuss 
the dynamics of charged dendrimers in a salt-free solu- 
tion. 
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FIG. 6: Radial number densities p(r) as a function of the distance r from the center of mass (CM) of the dendrimer for Ab 
equal to 0.0; 0.2Z ; 0.8 1 ; and 8.0 Z . 




FIG. 7: The number of beads (corresponding to the number densities in Fig. |SJ as a function of the distance r from the CM 
of the dendrimer. 



A. Global reorientational dynamics 



First, we consider the relaxation of vectors pointing 
from the center of mass (CM) of the dendrimer to the 
beads of a given generation shellii2i2L2£ If the terminal 
(G4) generation shell is considered, the corresponding 
autocorrelation function (ACF) of "CM-G4" vector can 
be employed to analyze the global reorientational behav- 



ior of the whole dendrimer^i It is, however, necessary 
to bear in mind that there is no reason for the beads 
of the last generation to be located at the dendrimer's 
periphery. 

In general, relaxation of the "CM-G4" vector is deter- 
mined by the rotation of the whole dendrimer together 
with fluctuations in the vector's length. To character- 
ize the dendrimer's rotation, we calculated the ACF of 
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FIG. 8: Radial charge densities p e (r) as a function of the 
radial distance r from the CM of a dendrimer for Ab equal 
to 0.2 lo (solid line), 0.4 lo (dashed line), 0.8 lo (solid line with 
circles), and 3.2 lo (solid line with stars). To reduce the noise 
in the data, the charge densities shown here were first fitted 



to splines.' 
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FIG. 9: Characteristic relaxation times r s h e ii and r ro t related 
to the dendrimer's reorientational dynamics as extracted from 
the corresponding autocorrelation functions of vectors point- 
ing from the CM of the dendrimer to the beads of its terminal 
G4 generation shell, see text for details. Error bars are of the 
order of the symbol size. 



the unit vector directed along the "CM-G4" vector. It 
turns out that the corresponding relaxation times r rot 
are very close to the relaxation time of the "CM-G4" 
vector (r s hoii)- Thus, the relaxation of the "CM-G4" vec- 
tor is governed almost exclusively by the rotation of the 
dendrimer in agreement with the results from Brownian 
dynamics simulations of dendrimers with rigid bonds^ 
Figure shows that the global reorientational motion 
of the charged dendrimer depends non-monotonically on 
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FIG. 10: Time evolution of the normalized autocorrelation 
function Pi(t) = (bj(O)bj(t)) of bond vectors bi(t) of the in- 
nermost (top) and the outermost (bottom) shells. The legend 
shows the different Bjerrum lengths used. 



the Bjerrum length. Such a behavior resembles the Ab- 
dependence of (Rg) as was seen earlier in Fig. El This is 
expected since it is reasonable to assume that the global 
reorientational dynamics is controlled by the dendrimer's 
size: larger size requires more time for reorientation. 

As Fig. |5| shows, there is a sharp maximum at Ab — 
0.8 lo. It is also interesting that although at large Bjer- 
rum lengths the average size of the charged dendrimer is 
found to be smaller than the size of a neutral dendrimer 
(Fig. |5J), that is not reflected in the characteristic relax- 
ation times shown in Fig. [5] The former can be explained 
by the inhomogeneous distribution of counterions inside 
the dendrimer, giving rise to the additional electrostatic 
attraction between different parts of the dendrimer in- 
terior. On the other hand, the global reorientation dy- 
namics is not only governed by the dendrimer's size but 
also by its mass. At Ab = 8.0 lo almost all counterions 
are bound to the dendrimer (Fig. 0}, thus increasing the 
dendrimer's effective mass by a factor of 1.5 in compar- 
ison to the neutral one. This, in turn, slows down the 
global reorientational motion of a dendrimer. 
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FIG. 11: The characteristic relaxation time r coro for a vector 
pointing from the CM of a dendrimer to its core as a function 
of the Bjerrum length Ab. Error bars are of the order of 
symbol size. 

B. Local dynamics of bonds and dendrimer's core 

As seen in Fig. |5J the global reorientational dynam- 
ics is determined by rather long characteristic times up 
to ^200 ro. The local orientational motions are charac- 
terized by much shorter times. As an example, we con- 
sider here the re-orientation dynamics of the bond vectors 
hi (t) . In Fig. EI] we plot the time evolution of the cor- 
responding autocorrelation function Pi(t) = (bj(O)bj(f)) 
for the innermost (GO) and the outermost (G4) shells. 
We note that the relaxation of Pi(t) cannot be charac- 
terized by a single relaxation time. Therefore, we focus 
here on the time behavior of P± (t) rather than on char- 
acteristic times. 

The reorientation of a bond slows down with the size 
of the dendrimer, see Fig. ^| (top). Such a behavior 
is found to hold for all inner generation shells from Gl 
to G3 (data not shown). Remarkably, the above non- 
monotonic (notice the order of the curves) dependence 
on Ab in Fig. E3 (top) breaks down for P\{t) of the ter- 
minal (outermost) generation shell as it is seen in Fig. 1101 
(bottom). This feature can be linked to condensation of 
counterions which plays a dominant role at Ab > 0.8 Iq 
and, therefore, can hinder the reorientation of bonds of 
the terminal G4 shell. 

The local re-orientation dynamics of dendrimer bonds 
can also be characterized by the autocorrelation function 
P 2 (t) = (3/2)((bi(0)b 4 (i)) 2 - 1/3). It turns out that 
the simple relation Piif) = P\(t) holds for all bonds in 
agreement with the previous simulation studies^ 

Last, we focus on the relaxation of the vector pointing 
from the dendrimer's center of mass to its core ("CM- 
core" vector) i& The core and the CM are typically found 
to be very close to each other. However, they do not 
exactly coincide and the relaxation of the corresponding 



"CM-core" vector is characterized by small (but finite) 
relaxation times. 

In Fig. II II we plot the corresponding characteristic re- 
laxation times r corc . As argued above, the relaxation 
times turn out to be very small, just a few tq. The Ab- 
dependence of the r corc resembles the behavior of the 
(-Rg)-curve in Fig. [3 i.e., increase in the dendrimer size 
speeds up the relaxation of the "CM-core" vector and 
vice versa. This effect can be explained in terms of the 
mobility of a dendrimer core: larger dendrimcrs are char- 
acterized by smaller densities of their interior (see Fig. [fjj) 
and, therefore, by more mobile cores which relax faster. 



V. DISCUSSION AND CONCLUSIONS 

We have performed a systematic molecular dynamics 
study of generic cationic dendrimers solvated in explicit 
solvent in the presence of explicit counterions. The main 
goal was to understand to what degree the structural and 
dynamic properties of charged dendrimers are sensitive to 
the strength of electrostatic interactions, and what is the 
role of counterions. 

We found an intriguing interplay between the repulsive 
interactions of the dendrimer's charged terminal beads 
and their attractive interactions with oppositely charged 
counterions. Depending on the type of dominant interac- 
tions, one observes swelling or shrinking of the molecule 
as the strength of the electrostatic interactions is in- 
creased. In other words, the size of a charged dendrimer 
depends non-monotonically on the Bjerrum length Ab 
as shown in Fig. [5] Most of the dendrimer's structural 
and dynamic characteristics considered in this study were 
found to follow this pattern. 

For small Bjerrum lengths a charged dendrimer swells 
with Ab and its size reaches maximum at Ab = 0.8 ^o- 
At larger Bjerrum lengths the condensation of counte- 
rions dominates and the dendrimer shrinks. The coun- 
terion condensation is often discussed in the framework 
of the Manning theory^ which predicts the condensa- 
tion to occur when Ab/o > 1 where a is the distance 
between charges along the polymer backbone. However, 
the Manning's criterion is not applicable in our case. The 
Manning's theory is formulated for linear polyelectrolytes 
and not for spherical-like branched objects such as den- 
drimers. Second, the charged dendrimer considered in 
our study has only its terminal beads charged, making 
the meaning of the parameter a uncertain. The counte- 
rion condensation, the effects of valence, charging, and 
the size and shape of the dendrimer could more appro- 
priately be studied using the method of Patra et al4£ 
who studied the charging of spherical three dimensional 
objects under different conditions. A detailed study of 
the above factors is, however, beyond the current work 
and will be addressed in a future publication. 

The Bjerrum length Ab = <J corresponds to the po- 
sition of the (iig)-peak, see Fig. |2 and, therefore, sep- 
arates swelling and shrinking regimes. Remarkably, on 
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the left-hand side of the peak, the separation between 
charges always exceeds the Bjerrum length because the 
beads cannot get closer than a and Ab < er. On the other 
hand, the counterion condensation is rather weak in this 
Ae-domain. Therefore, one can expect the Debye-Huckcl 
approximation to be valid at around Ab < cr and to break 
down at larger Bjerrum lengths. 

It is instructive to relate the dendrimer bond length Iq 
(and, correspondingly, the average dendrimer size i? g and 
the Bjerrum length Ab) to the characteristic length scales 
typical for real dendrimers in aqueous solutions. The 
value of 17.1 A is available from small angle X-ray scat- 
tering experiments^ for the average radius of gyration 
R g of a PAMAM dendrimer of the fourth generation. A 
very recent MD study2I reported R g = 16.78±0.15 A for 
a non-protonated PAMAM dendrimer in explicit water. 
In another study, Lee et aliSi reported R g = 14.8 ±0.1 A. 
That was obtained from MD studies of PAMAM den- 
drimers without solvent (R g = 14.50 ± 0.28 A in Ref. l2ll 
was obtained under the same conditions). Assuming that 
dendrimer swells on about 15 % in water the results of 
Lee et al~ seem to be in accord with other studies. We 
took R g ~ 17 A for our parameterization. For a neutral 
dendrimer we have i? g ~ 2.54 Iq, see Fig. [21 such that 
l ^ 6.7 A. 

Thus, the critical Bjerrum length corresponding to the 
dendrimer 's maximal size can be estimated to be 5.4 A. In 
turn, the Bjerrum length in water at 300 K is about 7.1 A 
(or ~ 1.3/o), i-e., it exceeds the critical value. In prac- 
tice, since the (i?g)-peak is not very sharp, the value of 
{Rg) for a dendrimer in water, being located next to the 
peak, is rather close to the maximum of {Rg), see Fig. [21 
Taking also into account the uncertainty in defining the 



value of Iq, one can conclude that counterion condensa- 
tion will not affect significantly PAMAM dendrimers in 
water under physiological conditions. 

To summarize, we demonstrated that including explicit 
counterions can have a dramatic effect on the structure 
and dynamics of charged dendrimers, and under certain 
conditions they cannot be treated implicitly. Therefore, 
simplified approaches, such as the well-known Debye- 
Hiickel approximation, have to be applied with great 
care. Given that the most exciting applications of proto- 
nated dendrimers are related to drug and gene delivery, 
it would be interesting to clarify the role of counterions 
in complexation of charged dendrimers with linear poly- 
electrolytes bearing opposite charges ! 49 ' 50 This problem 
is addressed in our ongoing studies. 
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